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ABSTRACT 

We present a determination of the distributions of gamma-ray photon flux - the so 
called LogiV-LogS' relation - and photon spectral index for blazars, based on the third 
extragalactic source catalog of the Fermi Gamma-ray Space Telescope’s Large Area 
Telescope, and considering the photon energy range from 100 MeV to 100 GeV. The 
dataset consists of the 774 blazars in the so-called “Clean” sample detected with a 
greater than approximately seven sigma detection threshold and located above ±20° 
Galactic latitude. We use non-parametric methods verified in previous works to recon¬ 
struct the intrinsic distributions from the observed ones which account for the data 
truncations introduced by observational bias and includes the effects of the possible 
correlation between the flux and photon index. The intrinsic flux distribution can be 
represented by a broken power law with a high flux power-law index of -2.43±0.08 and 
a low flux power-law index of -1.87±0.10. The intrinsic photon index distribution can 
be represented by a Gaussian with mean of 2.62±0.05 and width of 0.17±0.02. We also 
report the intrinsic distributions for the sub-populations of BL Lac and FSRQ type 
blazars separately and these differ substantially. We then estimate the contribution 
of FSRQs and BL Lacs to the diffuse extragalactic gamma-ray background radiation. 
Under the simplistic assumption that the flux distributions probed in this analysis 
continue to arbitrary low flux, we calculate that the best fit contribution of FSRQs is 
35% and BL Lacs 17% of the total gamma-ray output of the Universe in this energy 
range. 

Key words: methods: data analysis - galaxies: active - galaxies: jets - BL Lacertae 
objects: general 


In this work we investigate the distributions of gamma-ray 
photon flux and photon index for the largest flux-limited 
gamma-ray sample of blazars available at the moment, that 
of the third extragalactic source catalog of the Fermi Larg e 
Area Telescope (LAT) (3LAC - lAckermann et al.l 1201^ 1. 
Most of the extragalactic objects observed by the Large 
LAT on the Fermi Gamma-ray Space Telescope are clas¬ 
sified as blazars (e.g. lAbdo et al.|l2010al ). the active galactic 
nuclei (AGNs) in which t he jet is pointing toward us (e.g. 
iBlandford &: Konig]||l979l ). Gharacterizing the intensity and 
spectral characteristics of the gamma-ray emission seen in 
blazars is an essential for understanding of the physics of the 


1 INTRODUCTION 


accre tion disk and black hole systems in AGNs (e.g. IPermerl 
I2OO7I) . Understanding the characteristics of blazars is also 
crucial for evaluating their contribution to the extragalactic 
gamma-ray background (EGB) radiation. 


* E-mail: jsingaiarichinond.edu 
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The Fermi-LAT source catalogs are based on observa¬ 
tions extending down to gamma-ray energies of 100 MeV. 
Although the 3LAC reports the photon fluxes in the range 
from 1 to 100 GeV, the photon flux in the range from 100 
MeV to 100 GeV is an important quantity for blazars be¬ 
cause the most photons are at the lower energies of this 
range, as are the most photons of the EGB. Therefore, we 
base this analysis on the full 100 MeV to 100 GeV photon 
flux, reconstructed as described in ^ If the whole 100 MeV 
to 100 GeV information is used, the threshold flux for detec¬ 
tion depends strongly on an object’s gamma-ray spectrum. 
This arises from the energy dependence of the point spread 
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function llAtwood et al.l[2009l i and is such that harder spec¬ 
tra are detected at lower fluxes. Thus, for determination of 
the flux distribution one needs to take into account both the 
flux and the photon spectral index, and that one then deals 
with a bi-variate distribution of fluxes and indexes, which 
is truncated because of the above mentioned observational 
bias. Because the flux detection threshold depends on the 
photon index, the observed raw distributions do not provide 
the true LogA^-LogS counts or the true distribution of the 
photon index. Thus a bias free determination of the distribu¬ 
tions is more complicated than just counting sources. Here 
we use non-parametric methods to determine the distribu¬ 
tions directly from the data at hand. Ot her methods such a s 
monte carlo have been used, such as in lAbdo et al.l (l2010cl l 
- hereafter MA, and, for example, with redshift information 
and suitable assumptions the distribution of photon index 
can be corrected for the b ias towards hard spectrum sources 
as in IVenters et al.l l|2009l). _ 

In a previous work (ISingal et al.l (l2012l l - hereafter BPl) 
we demonstrated the techniques with real and simulated 
blazar flux and photon index data from the first year Fermi- 
LAT extragalactic source catalog (ILAC - lAbdo et al.l 
l2010bl '). As explored by e.g. Petrosian (1992), when dealing 
with a bi-(or more generally multi)-variate distribution, one 
must determine the correlation (or statistical dependence) 
between the variables, which cannot be done by simple pro¬ 
cedures when the data are truncated. We use a method 
based on t he techniques first developed b y Efron and Pet¬ 
rosian (EP: lEfron fc Petrosianiri992l . ll999l l which determine 
the intrinsic correlations (if any) between the variables and 
then the mono-variate distribution of each variable, account¬ 
ing for the truncations of the data. These techniques have 
been proven useful for application to many sources with var¬ 
ied characteristics and most recently to the gamma-ray flux 
and photon index of blazars i n BPl and to radi o and op¬ 
tical luminosity in quasars in ISingal et al.l (120131 '). where a 
more thorough discussion and references to earlier works are 
presented. 

In ISingal et al.l (l2014l l we used the spectroscopic red- 
shift information provided for the complete sample of ILAC 
Flat Spectrum Radio Quasar (FSRQ) type blazars to deter¬ 
mine the gamma-ray luminosity function, its evolution, and 
the density evolution for FSRQs. A very similar sample was 
used with different analysis techniques in i Aiello et al.[ (l2012l l 
to obtain the gamma-ray luminosity function and evolution 
for FSRQs, and a ILAC-based sample with a combination 
of measured redshifts and upper and lower limits w as used 
to do the same with BL Lacs in lAiello et al.l (l2014l L How¬ 
ever, in the case of 3LAC blazars, the spectroscopic redshift 
information is incomplete, leading to the utility of this work 
recovering the intrinsic flux and photon index distributions 
using the maximal number of well characterized gamma-ray 
detected blazars of any type. 

In this paper we apply these methods to determine the 
correlation and the intrinsic distributions of flux and pho¬ 
ton index of Fermi-LAT blazars. In (J^we discuss the data 
used from the 3LAC extragalactic catalog, and in !j3]we ex¬ 
plain the techniques used and present the results. In ^we 
describe how the results from such studies are important 
for understanding the origin of the extragalactic gamma-ray 
background (EGB) radiation. A discussion is presented in 

§ 5 . 
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Figure 1. Flux and reported photon spectral index for the 774 
Fermi-LAT blazars used in this analysis, those with test statistic 
^ 50 and |fe| IS 20°. BL Lac type blazars (n=352) are shown as 
blue triangles, FSRQs type blazars (n=286) are shown as red plus 
signs, and blazars of unidentified or ambiguous type (n=136) are 
represented by black x’s. It is seen that there is a selection bias 
against soft spectrum sources at fluxes below ~ 6 X 10“® photons 
cm~^ sec“^. We also show for a selection of sources (but only a 
few for clarity) the approximate limiting flux for that source - 
that is the lowest flux it could have and still be sufficiently bright 
to be included in the sample given its location on the sky given 
the reported detection significance. The location of the line used 
for the truncation boundary is also shown. 

2 DATA 

For this analysis we use the blazar sources reported in the 
Fermi-LAT 3LAC l|Ackermann et al]l2015al ) that are part of 
the “Clean” sample, have a detection test statistic TS ^ 50 
and which lie at Galactic latitude |6| ^ 20°. This set of cri¬ 
teria, which has been adopted by the LAT team for certain 
analyses of the blazar population (e.g. MA), includes those 
sources that are fully calibrated, removes spurious sources, 
and maximizes the likelihood of sources being properly iden¬ 
tified as to spectral type. The test statistic is defined as 
TS = —2 X (In(Lo) — In(Li)), where Lq and Li are the 
likelihoods of the background (null hypothesis) and the hy¬ 
pothesis being tested (e.g., source plus background). The 
significance of a detection is approximately n x a = y/TS . 

These include 286 which are identified as FSRQ type, 
352 which are identified as BL Lacertae (BL Lac) type, and 
136 which have uncertain type. This can be compared to 
the ILAC with 352 TS ^ 50 and |fe| ^ 20° blazars, of which 
161 are FSRQs, 163 are BL Lacs, and 28 which are identi¬ 
fied as uncertain type. Thus the dataset used in this analy¬ 
sis contains more than twice as many blazars as that used 
in the previous analyses of BPl and MA. There have also 
been some blazars that have been reclassified as a different 
ty pe since the ILAC , base d on e.g. sp e ctral data obtained 
by IShaw et al.l ll2012l) and IShaw et al.l ll2013l) among other 
works. 

The gamma-ray fluxesof the blazars in the range from 
100 MeV to 100 GeV, designated here as Fioo, and reported 
photon indexes are plotted in Figure [T] The 774 blazars 
range in Fioo from 6.01 x 10“^° to 2.10 x 10“® photons 
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cm“^ sec“^. The photon index T is defined such that for the 
photon spectral density at energy E, n{E)dE oc E~^ (or the 
uEi, oc and is obtained by the LAT collaboration 

by fitting a power-law to the observed spectra in the above 
energy interval. In this dataset the reported photon index 
ranges from 1.438 to 3.100. It should be noted that some 
blazars manifest a photon flux as a function of photon en¬ 
ergy in the range from 100 MeV to 100 GeV that cannot be 
completely described by a single power law, showing curva¬ 
ture in the spectrum or other features. However, the 3LAC 
does report a best-fit power law index for every blazar, and 
this analysis uses that power law. 

We recover Tioo from the reported flux density (K), 
pivot energy (Ep), and photon index with 

/•lOOGeV / \ r 

Fioo= K[—]dE (1) 

J lOOMeV \^P / 

Flop is the blazar g amma-ray flux measure adopted in e.g. 
lAbdo et al.l (|2010d l and is useful for exploring blazars and 
their relationship to the EGB as discussed in m We obtain 
the flux density and the pivot energy for 3LAG blazars from 
the c orresponding general Fermi source catalog dAcero et al.l 
l2015h . The bias mentioned above is clearly evident; there is 
a strong selection against soft spectrum sources at fluxes 
below Flop ~ 6 X 10~® photons cm“^ sec“^. 

Every blazar has an associated TS as discussed above, 
with the background flux b eing a function of position on 
the sky (lAbdo et al.l I2OIOIJI . In Figure [T] we also show the 
approximate limiting flux of some (only some to avoid vi¬ 
sual confusion) of the blazars, which is an estimate of the 
lowest flux it coul d have t o be included in the sample, given 
by Fum ~ Fipp/y^TS/50. This Fum will thus have some de¬ 
pendence on the object’s photon index and position on the 
sky. However, as discussed in BPl, since there is some uncer¬ 
tainty in the detection threshold values of fluxes and indexes, 
the limiting flux as determined in this way is not the optimal 
estimate. We therefore use a more conservative truncation 
as shown by the straight line in Figure [T] As verified with 
the simulated data discussed in BPl, moving the truncation 
line to the right and down eliminates more sources in the 
edges of the sample, and does not change the results but in¬ 
creases the uncertainty, while moving it to the left changes 
the results. We choose the truncation line boundary that is 
at the edge of this region where the results are not changed. 
This way we lose some data points but make certain that 
we are dealing with a complete sample with a well defined 
truncation. 


3 DETERMINATIONS OF DISTRIBUTIONS 
3.1 Correlations 

As discussed in |JT] with a bi-variate truncated dataset one 
must determine whether the variables are independent be¬ 
fore considering distributions. If in this case Fioo and F are 
independent, the combined distribution 6 (^ 100 , F) can be 
separated into two independent distributions iH^Fipp) and 
h(r). However, they may not be independent even though 
flux clearly depends strongly on redshift while photon in¬ 
dex may not, and in fact does not in the case of FSRQs 
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Figure 2. Test statistic r versus intrinsic correlation factor /3 for 
a photon index and flux correlation of the form given in Equation 
[S] for the 3LAC blazars (solid curve), the subset of BL Lac type 
blazars (dashed curve), and the subset of FSRQs type blazars 
(dash-dot curve). The Itr range of best fit values for /3 are where 
|t| ^ 1. For comparison, the dotted curve shows the correlation 
factor for just those sources above 6 X 10~® photons cm“^ sec~^, 
where the data truncation in the Fioo,r plane is not as relevant. 


dSingal et al.l l2014l) . Correlations between flux and pho- 
ton index are observed for examp l e in gamma-ray bursts 
ijYonetoku et al.l l2004l : iLlovd et al.l l200flh . There may also 
be an observed correlation induced by the selection biases, 
as there certainly is in this case as can be seen in Figure [T] 
which should be removed in order to obtain bias free distri¬ 
butions of the variables. Thus, the first task is to establish 
whether the variables are independent. 

While the flux and photon index exhibit a strong corre¬ 
lation in the raw (and heavily biased) data, determining the 
intrinsic correlation when the data are truncated requires 
statistical methods to account for the missing data. In BPl 
we discuss in detail the methods we apply to the bi-variate 
Fermi-LAT observational blazar data. We use a technique 
first e xplor ed by Efron and Petrosian (lEfron fc PetrosianI 
1 1991 1 19991) to determine whether the two variables are 
correlated. This method utilizes a modified version of the 
Kendall Tau test with the test statistic 




( 2 ) 


to quantify the independence of two variables in a dataset, 
say {xj,yj) for j = 1,... ,n. Here TZj is the dependent vari¬ 
able (y) rank of the data point j in a set associated with it. 
For untruncated data (i.e. data truncated parallel to the 
axes) the set associated with point j includes all of the 
points with a lower (or higher) independent variable value 
{xk < Xj). If the data is truncated the unbiased set is then 
the associated set consisting only of those points of lower(or 
higher) independent variable (x) value that would have been 
observed if they were at the x value of point j given the trun¬ 
cation. 

If {xj,yj) are uncorrelated then the ranks of all of the 
points TZj in the dependent variable within their associ- 
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ated set should be distributed uniformly between 0 and 
the number of points n, with the rank uncorrelated with 
their independent variable value, with the expectation value 
Sj = (l/2)(n + 1) and variance and Vj = (l/12)(n^ — 1), 
where n is the number of objects in object j’s associated 
set. Then the points’ contributions to r will tend to sum 
to zero. On the other hand, if the indepenent and depen¬ 
dent variables are correlated, then the rank of a point in the 
dependent variable will be correlated with its independent 
variable, and because the set to be ranked against consists 
of points with a lower independent variable value, the con¬ 
tributions to r will not sum to zero. 

Independence of the variables is rejected at the m a level 
if I T I > m, and this can be considered the same standard 
deviation as would be calculated from an other method such 
as lea st-squares fitting, as discussed in lEfron fc Petrosi^ 
(Il999ll . If the variables are not independent, to find the best 
fit correlation the y data are then adjusted by dehning some 
new dependent variable y'j = F{yj,Xj) and the rank test is 
repeated, with different values of parameters of the function 
F. 

To determine the intrinsic correlation between fioo and 
r we use a function which is a simple coordinate rotation, 
defining a new variable we call the “correlation reduced pho¬ 
ton index” as 

r„ = r-^ X log(:^) . (3) 

Then we determine the value of the parameter /I empiri¬ 
cally that makes Eioo and Ter independent. This is the best 
fit value of the correlation between Tioo and T for this func¬ 
tional form. The distributions of Fioo and Ter are indeed 
separable: 

G{Fioo,r) = iiFioo) X U.r,r). (4) 

Once the monovariate distributions of Tioo and Ter are 
determined then the true distribution of T can be recovered 
by an integration over Tioo as: 


m = 


' Aoo 


V<i^ioo) h (r - /? X log (^)) rf^’ioo. 


(5) 


Here fb is some fiducial flux we chose to be To = 6 x 10““ 
photons s“^ cm“^ sr“^ which is approximately where the 
flux distribution breaks (see below), although its actual 
value is not important - it is simply utilized to avoid tak¬ 
ing the logarithm of a dimensionless number. The data de¬ 
scribed in 31] are truncated in the fioo-r plane, due to the 
bias against low flux, soft spectrum sources. As discussed 
in 31] we can use a curve approximating the truncation, 
riim = 9(logTioo), which defines the associated set for each 
point as those objects whose photon index is less than the 
limiting photon index of the object in question given the 
truncation curve and its specific value of Tioo. 

As discussed in BPl we have tested this procedure using 
a simulated dataset from the Temii-LAT collaboration that 
resemble the real observations and are subject to a trunca¬ 
tion similar to the actual observational data, but with known 
input distributions of uncorrelated photon index and flux. 
Note that when defining new variables the truncation curve 


as a function of flux should also be transformed by same 
parameter /I; 

rcr.iim = Tita -p X Log { ) . (6) 

—min J 

Figure [2] shows the value of the test statistic r as a func¬ 
tion of the correlation parameter P for a all blazars and the 
subsets including only BL Lacs and FSRQs in the sample. 
Table 1 shows the best fit values and la ranges of the corre¬ 
lation parameter /3 for these cases. We note that the correla¬ 
tion is generally small, for example for all Fermi blazars the 
best fit value is (?=-0.13 ± 0.06, but significant. This is not 
surprising, as there is some evidence of an invers e correlation 
betw een blazar luminosity and photon index ll Aiello et alJ 
l2014h . In ISingal et al.l (l2014l j we found no evidence of signif¬ 
icant evolution of photon index with redshift for the case of 
FSRQs, so the result of a possible negative correlation seen 
here between photon index and flux would indicate that the 
(small) intrinsic correlation is actually indeed between pho¬ 
ton index and luminosity. 


3.2 Distributions 

With the correlation removed the independent distributions 
i^Fioo) and h(Tci) can b e determined using a method out- 
lin ed in j Petrosi anl (1 1992] ) based on techniques first devised 
by iLvnden-Belll ( igTll L These methods give the cumulative 
distributions by summing the contribution fro m each point 
witho ut binning the data, and as discussed in iLvnden-Belll 
(I1971II . are equivalent to the maximum likelihood distribu¬ 
tions. 


3.2.1 flux distributions 

For the flux, the cumulative distribution 

poo 

^Fioo)= / V<^i'oo) rf^’i'oo (7) 

-P’100 

is obtained with 

<KR»)= n (1 + 74) 

where j runs over all objects with fluxes Tioo,j Tioo, 
and N{j) is the number of objects i in the associated set 
of object j that have a value of Fioo,i ^ Fioo,j . The asso¬ 
ciated set for object j in this case is all of those objects 
that have Fcr ^ rcr.iim(Tioo,j), determined from the trunca¬ 
tion curve described above. Equation [S] represents the estab¬ 
lished Lynden-Bell method which we modify by including 
only those objects within object j’s associated set for the 
calculation of N{j). The use of only the associated set for 
each object removes the biases introduced by the truncation 
in this calculation. 

We determine the differential flux distribution 

iiF.00) . (9) 

by fitting piecewise polynomial functions via least-squares 
fitting to $(Tioo) and calculating the derivatives. Figure 
[3] shows the calculated intrinsic differential distribution 
tflFioo), along with those obtained from the raw observed 
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Figure 3. Observed (diamonds) and reconstructed intrinsic 
(stars) differential distribution of gamma-ray flux i/^Fioo) for the 
3LAC Fermi-LAT blazars. The vertical error bars on the stars 
represent the Icr range of the correlation parameter /3, the dom¬ 
inant source of uncertainty, and are generally smaller than the 
points. The intrinsic distribution is a power law with a break at 
Fbr — 8 X 10“® photons cm~^ sec“^. The best fit slopes for 
the intrinsic distribution are -2.43ih0.08 above the break and - 
1.87ib0.10 below, and the best fit intrinsic distribution is plotted 
as the dotted line. We also plot the best-fitting iH^Fioo) as de¬ 
termined in BPl (small crosses) and MA with error bars (dotted 
lines). The best fit value for ip(F},f.eak) is 1.66 xlO^ sr“^_F^’^Q. 
The distribution parameters for the BL Lac and FSRQ sub¬ 
populations separately are presented in Table I. 



Figure 4. Observed (diamonds) and reconstructed intrinsic 
(stars) distribution of photon index h(r) for the 3LAC Fermi- 
LAT blazars used in this analysis. The intrinsic distribution is 
calculated from the flux distribution and the correlation reduced 
photon index distribution by equation [5] The stars represent the 
intrinsic distribution calculated with the best fit value of the cor¬ 
relation parameter /3 and the solid curve is the best fit Gaussian 
function to these values, while the dotted curves represent the 
best fit Gaussian functions to the extremal intrinsic distributions 
allowed by the la range of /3. The intrinsic distribution can be 
represented by a Gaussian with a mean of 2.62di0.05 and la width 
of 0.17di0.02, while the raw observed distribution can be repre¬ 
sented by a Gaussian with a mean of 2.19di0.01 and la width of 
0.34±0.01. The normalization of h(r) is arbitrary. 


data without correcting for the bias. A direct comparison 
to the results from BPl and MA is presented there as well. 

For consistency with earlier works such as MA we fit the 
differential counts as a broken power law which describes 
the data well and takes the form 

/ Finn \'’^o.bove 

HFlOo) = HFhreak) ( - 1 for FlQQ ^ Fbreak (10) 

N ^break / 

/ Finn \'^below 

iiFbreak) ( - 1 for -^100 < Fbreak- 

\ -C break / 


3.2.2 Photon index distributions 

To determine the cumulative distribution of the correlation 
reduced photon index 

CT 

P{r,r)= / h(Odr', (12) 

Jo 

we use 

= n (' + «fe) 


mabove and mbeiow are the power law slopes above and be¬ 
low the break, respectively, and are obtained from a least- 
squares fitting of iliFoo), as is the value of Fbreak- At values 
of Fioo above Fnt = 6 x 10~® photons cm“^ sec~^ the 
truncation imposed by the boundary line shown in Figure 
[T]is not significant and we can obtain the normalization by 
scaling the cumulative distribution $(-Fioo) such that 




n±Vn 

8.26 sr ’ 


( 11 ) 


where N is the number of objects above Fnt in the dataset 
in question and is equal to 89 for all blazars, 20 for BL Lacs, 
and 65 for FSRQs. The y/N uncertainty arises because of 
Poisson noise for the brightest sources, and 8.26 sr is the to¬ 
tal sky coverage considered, which is |6| ^ 20° as discussed 
in !j2] This normalization then allows the properly normal¬ 
ized value of 4{Fbreak) to be calculated. The best fit value for 
'4i.Fbreak), Corresponding to the best-fitting value of ruabove 
is then 1.2 xlO® sr“'^F’i7in. 


in the modified method of iLvnden-Belll ([TotJ) as above, 
which can be differentiated to give the differential distri¬ 
bution 


HTer) 


dP{Ter) 

dFcr 


(14) 


In this case, k runs over all objects with a value of rcr,k ^ 
Per, and Af(/c) is the number of objects i with P cr,i ^ Bcrjfe 
in the associated set of object k. Here the associated set for 
object k is those objects with Tioo ^ Tiim.k obtained from 
the truncation line at Fcr,k. We note again that the cutoff 
curve as a function of flux is scaled by /? in the same manner 
of equation [S] 

HVct) can be used via equation[S]to obtain the intrinsic 
distribution of the photon index itself, HT). The results are 
shown in Figured along with the raw observed distribution 
for comparison. Because the mean of intrinsic distribution of 
photon index is sensitive to the value of the correlation pa¬ 
rameter /3, we include the full range of intrinsic distributions 
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Fioo jphotons cm ® s 

Figure 5. Estimate of the cumulative number of photons be¬ 
tween 0.1 and 100 GeV above a given -Fioo? 2^y-pop(> Fioo) from 
equation 1161 from FSRQs (red upper points) and BL Lacs (blue 
lower points), shown with error bars resulting from the la range of 
the correlation parameter /3. The error bars are generally smaller 
than the plotting symbols used, except for the leftmost points. 
The dashed horizontal line shows the EGB, ie the total extra- 
galactic gamma-ray output (Tegb) defined here. An estimated 
total contribution of a population to Xegb can be obtained by 
integrating equation 1161 to zero flux, as discussed in aa 

resulting from the la range of /3. A Gaussian form provides 
a good description of the intrinsic distribution of the index. 

We have carried out the same procedures to obtain the 
distributions of the BL Lac and FSRQ subsets of the data. 
Table 1 summarizes the best fit parameters for the intrin¬ 
sic flux and photon index distributions, for the sample con¬ 
sidered as a whole, and for the BL Lac and FSRQs sub¬ 
populations separately. The errors reported include statis¬ 
tical uncertainties in the fits and the deviations resulting 
from the Icr range of the correlation parameter /3, with the 
later being much larger in magnitude. A higher value of /3 
(i.e. more positive correlation between flux and photon index 
absolute value) moves the mean of the photon index distri¬ 
bution down to a lower absolute value of the photon index 
and makes the faint end source counts slope less steep (less 
negative mheiow), while a lower value of /3 has the opposite 
effect. 


4 TOTAL OUTPUT FROM BLAZAR 
POPULATIONS AND THE 
EXTRAGALACTIC GAMMA-RAY 
BAGKGROUND 

One can integrate the flux distribution to calculate the to¬ 
tal flux from a population and its proportional contribution 
to a photon background. In this case, we are interested in 
the EGB, defined here as the total extragalactic gamma-ray 
photon outputQ In BPl we focused on using the flux distri- 

^ This definition avoids the problem that individual instruments 
resolve a different fraction of sources, leading to different esti¬ 
mates for the fraction of the total extragalactic photon output 
that is unresolved. 


bution to estimate the contribution from blazars as a whole. 
Here with a much larger sample we will use the flux dis¬ 
tributions to estimate the FSRQ and BL Lac contributions 
separately. The total output in gamma-ray photons from 
blazar sources with fluxes greater than a given flux Fioo, in 
terms of photons s“^ cm“^ sr“^ between 0.1 and 100 GeV, 
is 

/ CO 

Floo^Fioo)dFioo- (15) 

loo 

Integrating by parts the contribution to the EGB can be 
related directly to the cumulative distribution <I>(Fioo) which 
is the primary output of our procedure 

poo 

l-y{> Fioo) = Fioo^Fioo)+ / $(F(oo)rfF'i'oo- (16) 

'' ^100 

The advantage of using the latter equation is that it can give 
a step-by-step cumulative total contribution to the back¬ 
ground instead of using analytic fits to the differential or 
cumulative distributions obtained from binning the data. 
Figure [S] shows Fry{> Fioo) for FSRQs and BL Lacs at fliLxes 
probed by this analysis. Note that this includes the contri¬ 
bution from both detected blazars and some of those un¬ 
detected owing to the truncation in the Fioo,r plane but 
still probed by the analysis. Therefore, this calculated con¬ 
tribution can be more than the total contribution of blazars 
resolved by Fermi-LAT in the SLAG. 

In order to estimate the contribution of objects at fluxes 
not probed by the Fermi-LAT to the total EGB one can 
extrapolate the flux distribution we have obtained to lower 
fluxes. This extrapolation is more uncertain. We fit a power 
law to the faint end of the "I^Fioo) distributions so that we 
can extend the integration of equation [T^ to lower fluxes. 
Integrating to zero flux we find that FSRQs would produce 
a photon output of F-y{> Fioo = 0)=3.6 (-I-3.4/-0.9) xl0“® 
photons s“^ cm“^ sr“^ and BL Lacs would produce a photon 
output oi X.y{> Fioo = 0)=1.8 (-I-4.5/-0.5) xl0“® photons 
s“^ cm“^ sr~^ . These relatively large ranges are due to the 
uncertainty in the faint end cumulative source counts slope, 
ultimately owing to the range of the correlation parameter 
P, where the best fit value reported is for the middle of the 
la faint end slope of "I^Fioo). 

This is to be compared with the total observed EGB, 
consisting of unresolved emission plus resolved sources, of 
Xegb = 10.4 x 10~^ photons s ~^ cm“^ sr~^ reported by 
Fermi (lAckermann et ani2015bh . It should be noted that 
this represents a change from the previous Fermi estimate 
of Xegb ~ 14.4 ± 1.9 x 10“® photons s“^ cm“^ sr~^ 
llAbdo et al.ll2010dh . //the populations continue to have the 
fitted power law distribution to zero flux then FSRQs would 
produce 35(-l-35/-9)% of the observed EGB and BL Lacs 
would produce 17(-|-44/-12)%. 

These results are roughly c onsist ent with a number of 
other findings. In ISingal et al.l (l2014li we consider redshift 
evolution effects with a smaller sample with complete avail¬ 
able spectroscopic redshifts and find that for FSRQs the 
contribution to the EGB is 22(-|-10/-4)%, with the result 
as reported here scaled by a factor of 1.38 for the ratio 
of the higher previously reported Fermi EGB intensity to 
the ne w recently reported lower intensity, while lAiello et aP 
ll2012li report 30.1(4-2.5/-1.7)% for FSRQs with the result 
as reported here likewise scaled. For BL Lacs, lAiello et al.l 
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(|2014) estimate that they contribute between 14%-20% of 
the EGB with the same scaling]. However, comparing the 
results obtained here from a large sample but with only flux 
information to those considering smaller samples with red- 
shift evolution information can be enlightening, as discussed 
in US] 

We note that blazars should not be present to arbitrar¬ 
ily low flux so the integration here to zero flux is an overesti¬ 
mate in particular for the high ranges of the above estimates 
which are most sensitive to the low flux integration limit. 
The best fit values obtained here do not favor blazars being 
the s ole significant contributo r to the EGB. Several a uthors 
(e.g. IStecker &: Venter^ I2OIII : lAbazaiian et al.l I2OI2I ) have 
suggested that blazars could be the primary s ource of the 
EGB. The spectral index of the EGB of ~2.4 (lAbdo et al.l 
l2010dll is consistent with the mean photon index of the 
blazars as determined in BPl and in MA (although less so 
with the mean photon index de termined here). In a similar 
vein, IVenters fc Pavildoul il201lh have shown that the spec¬ 
trum of the EGB is consistent with a blazar origin. Other 
possible source populations include starforming galaxies, 
which have been rec ognized as a possib l e maj o r contributor 
to the EGB by e.g. Stecker fc VentersI ll2011Il . [Fields et al.l 
ll2010ll . andjLaeki_et_ay (l201l|) , although this has be en dis¬ 
puted by iMakiva et al.l (l201ll ~) , radi o galaxies (e.g Inou^ 
2011), and other non-blazar AGN fe.g. llnoue fc Totani 20091 . 


201111 . Additional arguments against blazars being the sole 


significant contributor to the EGB have been made b ased 
on the observed EGB anosotr opy in iGuoco et al.l (l2012ll and 
iHarding fc AbazadiiarJ (|2012| 1. 


5 DISCUSSION 

We have applied a method to calculate the intrinsic distri¬ 
butions in flux and photon index of blazars directly from 
the observed Fermi-LAT 3LAC catalog without reliance on 
simulations. This method accounts for the pronounced data 
truncation introduced by the selection biases inherent in the 
observations when the full 100 MeV to 100 GeV observa¬ 
tional range is used, and addresses the possible correlation 
between the variables. The accuracy of the methods used 
here when applied to Fermi-LAT blazar flux and photon in¬ 
dex data were demonstrated in the Appendix of BPl using 
a simulated dataset with known distributions. A summary 
of the best fit correlations between photon index and flux, 
and the best fit parameters describing the inherent distri¬ 
butions of flux and photon index, are presented in Table 
1 along with the values obtained in BPl with the smaller 
ILAG dataset and by MA. We have obtained the intrinsic 
distributions considering the major data truncation arising 
from Fermi-LAT observations. More subtle issues affecting 
the distributions we have derived, especially the photon in¬ 
dex distribution, may arise due to the finite bandwidth of 
the Fermi-LAT and lack of complete knowledge of the ob¬ 
jects’ spectra over a large energy range and deviations from 
simple power laws. However the Fermi-LAT bandwidth is 
sufficiently large that the contribution of sources which peak 


^ The original reported numbers are 16(- H0/-4)% for FSRQ s 
from BPl, 21.7(4-2.5/-1.7)% for FSRQs frorn lAi ello et al.l ll201 2h. 
and 10%-15% for BL Lacs from lAiello et all 1201^1 7^ 


outside of this range to the source counts and the EGB in 
this energy range will be small. 

A limitation of these results is the use of the reported 
best-fit power-law photon index for each blazar, even though 
some, or perhaps most, sources actually deviate from a strict 
power law in the rel ation between photon flux and photon 
energy. According to lAckermann et al.l ll2015ah a total of 91 
FSRQs, 32 BL Lacs, and 8 blazars of unknown type in the 
3LAC show significant curvature in their spectra. As shown 
in that paper it is overwhelmingly the highest TS (and there¬ 
fore generally higher flux) blazars which have been shown to 
have significant curvature away from a true power law in this 
relationship. Spectral curvature could be a property of the 
higher flux blazars, or possibly more likely it is a general 
feature of blazars which has been most readily observable in 
those blazars with the highest TS. If the former is the case, 
then the effect on the truncation and issues related to the 
truncation would be minimal no matter how the truncation 
is dealt with. However, even in the case of the later, since 
this analysis relies on an empirical curve to approximate the 
truncation in the reported Fioo,r plane as discussed in il3.1l 
the analysis itself should not be particularly sensitive to de¬ 
viations from a strict power law. However, to the extent that 
power laws are not representative of the true spectral prop¬ 
erties of bright blazars in the energy range from 100 MeV 
to 100 GeV, that will be the case in these results and the 
distributions of photon index must be understood to be dis- 
tributions o f the bes t fit power law photon index. As stated in 
lAbdo et al.l ll2010ell . “Although some spectra display signif¬ 
icant curvatures, the photon index obtained by fitting single 
power-paw models over the whole LAT energy range pro¬ 
vides a convenient means to study the spectral hardness.” 

In §3.2 of BPl we consider the uncertainties and errors 
present in this analysis resulting from individual measure¬ 
ment uncertainties, blazar variability, and source confusion. 
We conclude that the errors and uncertainties introduced 
from each of these effects are small compared to the dom¬ 
inant source of uncertainty in this analysis, which is that 
arising from the range of the correlation parameter fi which 
propagates through every other determination. As we argue 
in BPl, to the extent that faint blazars are more likely to 
be observed by instruments such as the Fermi-LPCL if they 
are in the flaring state ra ther than the quiescent state (e.g. 
IStecker fc Salamonlll996ll . then the observed blazars should 
have a different mean photon index than the EGB, were the 
EGB to be made primarily from quiescent state blazars, un¬ 
der the assumption that blazars in the flaring state have a 
different spectrum than in the quiescent state. As the re¬ 
constructed mean photon index here (and elsewhere) of the 
Fermi-LAT observed population is close to (although not 
exactly) that of the EGB, and there is only a weak relation 
and correlation between flux and photon index, this would 
imply that at least one of the following must be the case: 
a) there is not a significant bias in the Fermi-LAT toward 
detecting blazars in the flaring state, b) quiescent blazars 
do not form the bulk of the EGB, or c) flaring and quies¬ 
cent blazars have, en masse, roughly the same photon index 
distributions. 

We find that the photon index and flux show a slight 
negative correlation. This correlation is greater than Icr sig¬ 
nificance for all blazars, and is quite significant in the case 
of the FSRQ sub-population. This correlation likely results 
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from an underlying correlation between photon index and lu¬ 
minosity, as discussed in 1 13.11 Viewing the flux and photon 
index distributions, the comparison of the intrinsic and raw 
observed distributions show clearly the substantial effects of 
the observational bias, as can be seen in Figures [3] and [I] 
The intrinsic differential counts can be fitted adequately by 
a broken power law and the photon index appears to have 
an intrinsic Gaussian distribution. 

We also find that in general the values reported here 
are consistent with those reported in BPl and MA for the 
power law slopes of the flux distribution '!/(F’ioo) above the 
breaks and the distributions of photon index HT). The al¬ 
lowed range of the correlation parameter j3 here and in BPl 
allows for wider uncertainty in these values in some cases 
than in MA. However, as expected, the uncertainty is gen¬ 
erally lower here than in BPl, as this sample has roughly 
twice as many sources. We note that in this work we find sys¬ 
tematically steeper power law slopes of the flux distribution 
il^Fioo) below the breaks than in BPl and MA. The differ¬ 
ences with MA are significant although the differences with 
BPl are not significant at the la level due to the larger un¬ 
certainties. We also note that the best-fitting widths of the 
Gaussian fits to the photon index distributions are system¬ 
atically smaller here than in BPl and MA for the cases of 
all blazars and BL Lacs, but not for FSRQs, and that the 
mean photon index for all blazars found here is higher than 
in the other works, and this is entirely the result of the mean 
photon index for FSRQs being higher. 

Using the bias free flux distributions we calculated the 
integrated contribution of FSRQs and BL Lacs to the EGB 
as a function of flux. We obtain this directly from the cu¬ 
mulative flux distribution which is the main output of the 
methods used. Using the simplistic assumption that the dis¬ 
tributions continue unchanged to arbitrarily low flux, the 
best fit contribution of FSRQs and BL Lacs to the total 
extragalactic gamma-ray radiation in the range from 0.1 to 
100 GeV is estimated to be 35% and 17% respectively. The 
significant uncertainties reported here for the source count 
slopes and the estimated contribution of blazars to the EGB 
are ultimately due to the allowed range of the correlation 
parameter /], as discussed at length in BPl. These results 
for the proportional contribution to the EGB are roughly 
consis tent with a number of other findings. In ISingal et al.l 
ll20l4l we consider redshift evolution effects with a smaller 
sample with available redshifts and find that for FSRQs the 
contribution to the EGB is 22(-|-10/-4)%, while 
(20U) report 30.1(-|-2.5/-1.7)% for FSRQs and 
1 20141 estimate that BL Lacs contribute between 10%-15% 
of the EGB, with all of these reported here scaled to account 
for the subsequently reduced estimated level of the EGB as 
discussed in ^ However, we note that one would not nec¬ 
essarily expect exact agreement between the estimates of 
total photon output derived here extrapolating the flux dis¬ 
tributions to arbitrarily low flux and those derived consid¬ 
ering the full redshift evolutions, for two primary reasons. 
One is that as the populations clearly h ave both luminos¬ 
ity and density evolutio n in redshift ie.g. ISingal et al]l2014l : 
lAiello et al.ll20lll2014ll . one would not expect the flux dis¬ 
tribution to remain constant at lower fluxes. The other is 
that blazars likely do not continue to arbitrarily low flux, 
indeed in BPl we argue that a reasonable lower cutoff flux 
for blazars would be around ~ 10“^^ photons cm“^ sec“^. 


AielJo_et_^ 

Aiello et al. 


The fact that we do find reasonable agreement between the 
estimates here and those arrived at considering redshift evo¬ 
lutions indicates not only that analyses seem to be consis¬ 
tent, but also that there may be a “conspiracy” where lumi¬ 
nosity and density evolution effects, as well as an absolute 
luminosity cutoff for blazars, and any effects from contribu¬ 
tions missed here because some 17% of 2 LAC blazars are of 
unclassified type, sum in such a way as to render extrapo¬ 
lations of the distributions recovered here to arbitrarily low 
fluxes appropriate. It also points to the utility of comparing 
this analysis, which is based on the statistically larger 3LAC 
sample which is the largest extant sample of gamma-ray de¬ 
tected blazars, with those of smaller samples that contain 
complete redshift information. 
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Table 1 - Best fitting parameters for Fermi-ljAT blazar intrins ic distributions, as calculate d in this work, and two analyses with the 
smaller Fermi-LAT ILAC - ours (BPl - ISingal et al.|[2012l‘) and MA llAbdo et al.ll2010ct) 




n 

/3“ 

'^above 

^break 

^below 



Blazars^ (this work) 

774 

-0.13±0.06 

-2.43±0.08 

8.0±0.2 

-1.87±0.10 

2.62±0.05 

0.17±0.02 

BlazarsS (BPl) 


352 

0.02±0.08 

-2.37±0.13 

7.0±0.2 

-1.70±0.26 

2.41±0.13 

0.25±0.03 

BlazarsS (MA) 


352 

- 

-2.48±0.13 

7.39±1.01 

-1.57±0.09 

2.37±0.02 

0.28±0.01 

BL Lacs (this work) 

352 

-0.11±0.06 

-2.44±0.09 

6.0 ±0.2 

-2.00±0.15 

2.18±0.02 

0.18±0.03 

BL Lacs (BPl) 


163 

0.04±0.09 

-2.55±0.17 

6.5 ±0.5 

-1.61±0.27 

2.13±0.13 

0.24±0.02 

BL Lacs (MA) 


163 

- 

-2.74±0.30 

6.77±1.30 

-1.72±0.14 

2.18±0.02 

0.23±0.01 

FSRQs (this woi 

k) 

286 

-0.25±0.05 

-2.42±0.09 

9.2±0.1 

-1.55±0.12 

2.75±0.06 

0.19±0.02 

FSRQs (BPl) 


161 

-0.11±0.06 

-2.22±0.09 

5.1±2.0 

-1.62±0.46 

2.52±0.08 

0.17±0.02 

FSRQs (MA) 


161 

- 

-2.41±0.16 

6.12±1.30 

-0.70±0.30 

2.48±0.02 

0.18±0.01 


“The correlation between photon index P and Log flux Fioo- See Equation [J] and ^3.11 A higher value of (3 
(i.e. more positive correlation between flux and photon index absolute value) moves the mean of the photon 
index distribution down to lower photon index absolute value (lower /i) and makes the faint end source 
counts slope less steep (less negative fnieiow)^ while a lower value of 0 has the opposite effect. 

^The power law of the intrinsic flux distribution -i/^Fioo) at fluxes above the break in the distribution. See 
Equation 1101 All values reported for this work include the full range of results and their uncertainties when 
considering the Icr range of 0. 

“The flux at which the power law break in V^Fioo) occurs, in units of 10~® photons cm~^ sec“^. We present 
the value even though the precise location of the break is not important for the analysis in this work. 

^The power law of the intrinsic flux distribution -i/^Fioo) at fluxes below the break. See Equation 1101 
“The mean of the Gaussian fit to the intrinsic photon index distribution h(r). 

■^The Icr width of the Gaussian fit to the intrinsic photon index distribution HX). 

^Including all FSRQs, BL Lacs, and those of unidentified type. 
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